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DNS  OF  MULTIPLICITY  AND  STABILITY  OF  MIXED 
CONVECTION  IN  ROTATING  CURVED  DUCTS 


LIQIU  WANG  AND  TIANLIANG  YANG 
Department  of  Mechanical  Engineering 
The  University  of  Hong  Kong,  Hong  Kong 


Abstract.  A numerical  study  is  made  on  the  fully-developed  bifurcation 
structxire  and  stability  of  the  mixed  convection  in  rotating  curved  ducts 
with  the  emphasis  on  the  effect  of  buoyancy  force.  The  rotation  can  be 
positive  or  negative.  The  fluid  can  be  heated  or  cooled.  The  study  re- 
veals the  rich  solution  and  flow  structures  and  complicated  stability  fea- 
tures. One  symmetric  and  two  symmetric/asymmetric  solution  branches 
are  found  with  seventy-five  limit  points  and  fourteen  bifurcation  points. 
The  flows  on  these  branches  can  be  symmetric,  asymmetric,  2-cell  and  up 
to  14-cell  structures.  Dynamic  responses  of  the  multiple  solutions  to  finite 
random  disturbances  are  examined  by  the  direct  transient  computation.  It 
is  found  that  possible  physically  realizable  fully-developed  flows  evolve,  as 
the  variation  of  buoyancy  force,  from  a stable  steady  multi-cell  state  at  a 
large  buoyancy  force  of  cooling  to  the  co-existence  of  three  stable  steady 
multi-cell  states,  a temporal  periodic  oscillation  state,  the  co-existence  of 
periodic  oscillation  and  chaotic  oscillation,  a chaotic  temporal  oscillation,  a 
subharmonic-bifurcation-driven  asymmetric  oscillating  state,  and  a stable 
steady  2-cell  state  at  large  buoyancy  force  of  heating. 


1.  Introduction 

We  study  the  fully-developed  bifurcation-driven  multiplicity  and  dynamic 
responses  of  multiple  solutions  to  finite  random  disturbances  numerically  by 
the  finite- volume/Euler-Newton  continuation  and  the  direct  transient  com- 
putation for  the  mixed  convection  in  ducts  of  square  cross-section  with  the 
streamwise  curvature,  the  spanwise  rotation  in  either  positive  or  negative 
direction,  and  the  wall  heating/cooling  [Fig.  1 with  (R,Z,tp)  as  the  ra- 
dial, spanwise  and  streamwise  directions,  respectively].  A positive  rotation 
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gives  rises  to  a Coriolis  force  in  the  cross  plane  (i?Z-plane)  directed  along 
positive  R-direction  and  vice  versa.  Such  flows  and  transport  phenomena 
find  their  application  in  sedimentation  field-flow  fractionation,  aerosol  cen- 
trifuges, rotating  power  machinery,  rotating  heat  exchangers,  centrifugal 
material  processing  and  material  quality  control,  medical  and  chromato- 
graphic devices,  etc. 

Early  works  on  the  rotating  curved  duct  flows  were  constrained  to  two 
simplified  limiting  cases  with  strong  or  weak  rotations.  Ludwieg  (1951)  de- 
veloped a solution  based  on  a momentum  integral  method  for  the  isothermal 
flow  in  a square  duct  with  a strong  spanwise  rotation.  Miyazaki  (1971,  1973) 
examined  the  mixed  convection  in  a curved  circular/rectangular  duct  with 
spanwise  rotation  and  wall  heating  by  a finite  difference  method.  Because 
of  the  convergence  difficulties  with  the  iterative  method  used,  Miyazaki’s 
work  was  constrained  to  the  case  of  weak  curvature,  rotation  and  heating 
rate.  As  well,  all  the  works  employ  a steady  model  for  the  fully  developed 
laminar  flow  with  a positive  rotation  of  the  duct.  Since  the  solution  is  only 
for  the  asymptotic  cases,  the  secondary  flow  revealed  by  these  early  works 
consists  of  only  one  pair  of  counter-rotating  vortices  in  the  cross-plane.  The 
interaction  of  the  secondary  flow  with  the  pressure-driven  streamwise  flow 
shifts  the  location  of  the  maximum  streamwise  velocity  away  from  the  cen- 
ter of  the  duct  and  in  the  direction  of  the  secondary  velocity  in  the  middle 
of  the  duct. 

More  comprehensive  studies  have  been  made  in  recent  years  by  Wang  & 
Cheng  (1996a)  and  Daskopoulos  & LenhofF(1990)  for  a circular  tube,  Mats- 
son  & Alfredsson  (1990,  1994)  and  Guo  & Finlay  (1991)  for  a high-aspect- 
ratio  rectangular  duct,  and  Wang  & Cheng  (1995,  1996b,  1997,  2001),  Wang 
(1997a,  b,  1999),  Selmi  et  a/.(1994)  and  Selmi  & Nandakumar  (1999)  for 
the  square  and  rectangular  ducts  with  a low-aspect-ratio.  All  the  works  are 
for  the  steady  fully  developed  flows.  Wang  & Cheng  (1996a)  developed  an 
analytical  solution  for  rotating  curved  flow  with  effect  of  heating  or  cool- 
ing which  allows  to  analyze  the  solution  structure.  Detailed  flow  structures 
and  heat  transfer  characteristics  were  examined  numerically  by  Wang  & 
Cheng  (1996b)  and  Wang  (1997a,  b,  1999).  The  rotating  curved  flows  were 
visualized  using  smoke  injection  method  by  Wang  & Cheng  (1995,  1997, 
2001).  Daskopoulos  & Lenhoff  (1990)  made  the  first  bifurcation  study  nu- 
merically under  the  small  curvature  and  the  symmetry  condition  imposed 
along  the  tube  horizontal  central  plane.  Matsson  & Alfredsson  (1990)  pre- 
sented the  first  and  comprehensive  linear  stability  analysis.  Matsson  & Al- 
fredsson (1994)  reported  an  experimental  study,  by  hot-wire  measurements 
and  smoke  visualization,  of  the  effect  of  rotation  on  both  primary  and  sec- 
ondary instabilities.  Using  a linear  stability  theory  and  spectral  method, 
Guo  &;  Finlay  (1991)  examined  the  stability  of  streamwise  oriented  vortices 
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to  2D,  spanwise-periodic  disturbances  (Eckhaus  stability).  Detailed  bifurca- 
tion structure  and  linear  stability  of  solutions  was  determined  numerically 
by  Selini  et  al.  (1994)  and  Selmi  &;  Nandakumar  (1999)  without  imposing 
the  symmetric  boundary  conditions. 

It  is  the  relative  motion  between  bodies  that  determines  the  perfor- 
mances such  as  friction  and  heat  transfer  characteristics.  The  duct  rotation 
introduces  both  centrifugal  and  Coriolis  forces  in  the  momentum  equa- 
tion describing  the  relative  motion  of  fluids  with  respect  to  the  duct.  For 
isothermal  flows  of  a constant  property  fluid,  the  Coriolis  force  tends  to  pro- 
duce vorticity  while  the  centrifugal  force  is  purely  hydrostatic,  analogous  to 
the  Earth’s  gravitational  field  (Wang  2001).  When  a temperature-induced 
variation  of  fluid  density  occurs  for  non-isothermal  flows,  both  Coriolis  and 
centrifugal-type  buoyancy  forces  could  contribute  to  the  generation  of  the 
vorticity  (Wang  2001).  These  two  effects  of  rotation  either  enhance  or  coun- 
teract each  other  in  a nonlinear  manner  depending  on  the  direction  of  duct 
rotation,  the  direction  of  wall  heat  flux  and  the  flow  domain.  As  well,  the 
buoyancy  force  is  proportional  to  the  square  of  the  rotation  speed  while  the 
Coriolis  force  increases  proportionally  with  the  rotation  speed  itself  (Wang 
1997b).  Therefore,  the  effect  of  system  rotation  is  more  subtle  and  compli- 
cated and  yields  new,  richer  features  of  flow  and  heat  transfer  in  general, 
the  bifurcation  and  stability  in  particular,  for  non-isothermal  flows.  While 
some  of  such  new  features  are  revealed  by  our  recent  analytical  and  nu- 
merical works  (Wang  1997a,  b,  1999,  Wang  & Cheng  1996a,  b),  there  is 
no  known  study  on  the  bifurcation  and  stability  of  mixed  convection  in 
rotating  curved  ducts. 

The  present  work  is  a relatively  comprehensive  study  on  the  bifurca- 
tion structure  and  stability  of  multiple  solutions  for  the  laminar  mixed 
convection  in  a rotating  curved  duct  of  square  cross-section  (Fig.  1).  The 
governing  differential  equations  in  primitive  variables  are  solved  for  de- 
tailed bifurcation  structure  by  a finite-volume/Euler-Newton  continuation 
method  with  the  help  of  the  bifurcation  test  function,  the  branch  switching 
technique  and  the  parameterization  of  arc-length  or  local  variable.  Tran- 
sient calculation  is  made  to  examine  in  detail  the  response  of  every  solution 
family  to  finite  random  disturbances.  The  power  spectra  are  constructed 
by  the  Fourier  transformation  of  temporal  oscillation  solutions  to  confirm 
the  chaotic  flow.  We  restrict  ourself  to  the  hydrodynamically  and  thermally 
fully-developed  region  and  two-dimensional  disturbances.  So  far,  a detailed 
3D  numerical  computation  of  flow  bifurcation  and  stability  is  still  too  costly 
to  conduct.  A 2D  model  is  still  useful  for  a fundamental  understanding  of 
rotating  curved  duct  flows.  However,  our  assumption  of  fully  developed  flow 
limits  our  analysis  to  the  one  preserving  the  streamwise  symmetry.  There 
may  be  further  bifurcation  to  flows  that  breaks  this  symmetry  and  that 
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cannot  be  found  in  the  present  work. 

2.  Governing  Parameters  and  Numerical  Algorithm 

Consideration  is  given  to  a hydrodynamically  and  thermally  fully  developed 
laminar  flow  of  viscous  fluid  in  a square  duct  with  the  streamwise  curva- 
ture, the  spanwise  rotation,  and  the  wall  heating  or  cooling  at  a constant 
heat  flux  (Fig.l).  The  geometry  is  toroidal  and  hence  finite  pitch  effect  is 
not  considered.  The  rotation  can  be  positive  or  negative  at  a constant  an- 
gular velocity.  The  duct  is  streamwisely  and  peripherally  uniformly  heated 
or  cooled  with  a uniform  peripheral  temperature.  The  properties  of  the 
fluid,  with  the  exception  of  density,  are  taken  to  be  constant.  The  usual 
Boussinesq  approximation  is  used  to  deal  with  the  density  variation.  The 
gravitational  force  is  negligible  compared  with  the  centrifugal  and  Coriolis 
forces. 

Consider  a non-inertial  toroidal  coordinate  system  ( R , Z,  (f>)  fixed  to 
the  duct  rotating  with  a constant  angular  velocity  about  the  O' Z'  axis, 
as  shown  in  Fig.  1.  We  may  obtain  the  governing  differential  equations,  in 
the  form  of  primitive  variables,  governing  fully-developed  mixed  convection 
based  on  conservation  laws  of  mass,  momentum  and  energy.  The  bound- 
ary conditions  axe  non-slip  and  impermeable,  streamwise  uniform  wall  heat 
flux  and  peripherally  uniform  wall  temperature  at  any  streamwise  position. 
The  proper  scaling  quantities  for  non-dimensionalization  are  chosen  based 
on  our  previous  experience  (Wang  & Cheng  1996b).  The  formulation  of 
the  problem  is  on  full  flow  domain  without  imposing  symmetric  boundary 
conditions  to  perform  a thorough  numerical  simulation.  The  readers  are 
referred  to  Wang  & Cheng  (1996b)  for  the  details  of  mathematical  formu- 
lation of  the  problem. 

The  dimensionless  governing  equations  contain  five  dimensionless  gov- 
erning parameters:  one  geometrical  parameter  a (the  curvature  ratio  de- 
fined by  a/Rc  , the  ratio  of  duct  width/height  a over  the  radius  of  the 
curvature  Rc,  representing  the  degree  of  curvature),  one  thermophysical 
parameter  Pr  (the  Prandtl  number,  representing  the  ratio  of  momentum 
diffusion  rate  to  that  of  the  thermal  diffusion),  and  three  dynamical  pa- 
rameters Dk,  LI  and  L2  defined  in  Wang  & Cheng  (1996b).  The  pseudo 
Dean  number  Dk  is  the  ratio  of  the  square  root  of  the  product  of  inertial 
and  centrifugal  forces  to  the  viscous  force  and  characterizes  the  effect  of 
inertial  and  centrifugal  forces.  LI  represents  the  ratio  of  the  Coriolis  force 
over  the  centrifugal  force,  characterizing  the  relative  strength  of  Coriolis 
force  over  the  centrifugal  force.  L2  is  the  ratio  of  the  buoyancy  force  over 
the  centrifugal  force  and  represents  the  relative  strength  of  the  buoyancy 
force.  A positive  (negative)  value  of  LI  is  for  the  positive  (negative)  rota- 
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tion.  A positive  (negative)  value  of  L2  indicates  the  wall  heating  (cooling). 
In  the  present  work,  we  set  a = 0.2  (typically  used  in  cooling  systems  of 
rotor  drums  and  conductors  of  electrical  generators)  and  Pr  = 0.7  (a  typi- 
cal value  for  air)  to  study  the  effects  of  three  dynamical  parameters  on  the 
multiplicity  and  stability.  While  results  regarding  the  effects  of  Dk  and  LI 
are  also  available,  we  focus  on  the  effects  of  L2  at  Dk  — 300  and  LI  = 28 
in  the  present  paper  due  to  limited  space. 

The  governing  differential  equations  are  discretized  by  the  finite  vol- 
ume method  to  obtain  discretization  equations.  The  discretization  equa- 
tions are  solved  for  parameter-dependence  of  velocity,  pressure  and  tem- 
perature fields  by  the  Euler-Newton  continuation  method  with  the  solution 
branches  parameterized  by  L2,  the  arc-length  or  the  local  variable.  The 
starting  points  of  our  continuation  algorithms  are  the  three  solutions  at 
Dk  = 300,  LI  = 28  and  L2  = 0 from  our  study  of  the  effects  of  D and 
LI.  The  bifurcation  points  are  detected  by  the  test  function  developed  by 
Seydel  (1994).  The  branch  switching  is  made  by  a scheme  approximating 
the  difference  between  branches  proposed  by  Seydel  (1994).  The  dynamic 
responses  of  multiple  solutions  to  the  2D  finite  random  disturbances  are 
examined  by  the  direct  transient  computation.  The  readers  are  referred 
to  Wang  & Yang  (2001)  and  Yang  & Wang  (2000)  for  the  numerical  de- 
tails and  the  check  of  grid-dependence  and  accuracy.  The  computations  are 
carried  out  on  the  Super  Computer  SP2  of  The  University  of  Hong  Kong. 

3.  Results  and  Discussion 

3.1.  SOLUTION  STRUCTURE 

The  bifurcation  structure  is  shown  in  Fig.2  for  L2  values  from  —20  up 
to  70  at  cr  = 0.02,  Pr  = 0.7,  Dk  = 300  and  LI  = 28.  In  Fig.2,  the 
radial  velocity  component  u at  r = 0.9  and  z = 0.14  (where  the  flow 
is  sensitively  dependent  on  L2;  r = R/a  and  2 = Z/a ) is  used  as  the 
state  variable,  enabling  the  most  clear  visualization  of  all  solution  branches. 
Three  solution  branches,  labeled  by  AST,  AS 2 and  S3  respectively,  are 
found.  Here,  S stands  for  symmetric  solutions  with  respect  to  the  horizontal 
central  plane  z = 0,  and  AS  indicates  that  the  branch  has  both  symmetric 
and  asymmetric  solutions. 

Branch  AS1  has  sixty-nine  limit  points  labeled  by  AS11  to  AS1G9, 
eleven  bifurcation  points  connecting  its  sub-branches  denoted  by  AST451-1 
to  AST451-11,  two  bifurcation  points  connecting  itself  to  AS2  labeled  by 
AS1AS2~ 1 and  ASl'452-1,  and  one  bifurcation  point  connecting  itself  to  S3 
denoted  by  AS1S3.  Branch  AS2  has  four  limit  points  labeled  by  AS21  to 
AS23 4.  Branch  S3  is  a symmetric  solution  branch  and  has  two  limit  points 
S31  and  S32.  The  location  of  fourteen  bifurcation  points  and  seventy-five 
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limit  points  is  available  in  Yang  (2001).  To  visualize  the  details  of  branch 
connectivity  and  some  limit/bifurcation  points,  the  locally  enlarged  state 
diagrams  are  also  shown  in  Pig.  2.  As  Fig.  2 is  only  ID  projection  of  12400 
dimensional  solution  branches,  all  intersecting  points  except  fourteen  bifur- 
cation points  should  not  be  interpreted  as  connection  points  of  branches. 

For  a large  | L2  | value  (L2  < -14.5  or  L2  > 63.1),  the  buoyancy 
force  dominant  the  mixed  convection.  There  is  unique  flow  and  tempera- 
ture field  for  a specified  value  of  L2  in  these  two  ranges.  Figure  3 illustrates 
the  secondary  flow  patterns,  the  streamwise  velocity  profiles  and  temper- 
ature profiles  at  L2  = —17  and  L2  = 65,  respectively.  In  the  figure,  the 
stream  function,  streamwise  velocity  and  temperature  are  normalized  by 
their  corresponding  maximum  absolute  values  | ip  |m0a:  ? wmax  and  tmax.  A 
vortex  with  a positive  (negative)  value  of  the  secondary  flow  stream  func- 
tion indicates  a counter-clockwise  (clockwise)  circulation.  The  readers  are 
referred  to  Wang  & Cheng  (1996b)  for  a detailed  discussion  of  the  flow 
structures  shown  in  Fig.  3 in  general,  their  relations  with  physical  mecha- 
nisms and  driving  forces  and  their  effects  on  the  flow  resistance  and  heat 
transfer  in  particular. 

For  a L2  value  in  —14.5  < L2  < 63.1,  however,  we  can  have  multiple 
solutions.  Figure  4 shows  typical  secondary  flow  patterns  of  six  solutions 
(thirty-nine  solutions  in  total)  at  L2  — —11.7.  It  is  observed  that  the  nonlin- 
ear competition  of  driven  forces  leads  to  not  only  a rich  solution  structure 
but  also  complicated  flow  structures.  Therefore,  the  mixed  convection  in 
rotating  curved  ducts  is  much  more  complicated  than  that  available  in  the 
literature. 

3.2.  STABILITY  OF  MULTIPLE  SOLUTIONS 

Recognizing  that  there  is  no  study  on  dynamic  responses  of  multiple  solu- 
tions to  finite  random  disturbances  in  the  literature,  a relatively  compre- 
hensive transient  computation  is  made  to  examine  the  dynamic  behavior 
and  stability  of  typical  steady  solutions  with  respect  to  four  sets  of  finite 
random  disturbances  with  d = 4%,  10%,  15%  and  40%  respectively.  Here,  d 
is  the  maximum  percentage  of  disturbing  value  over  the  steady  value  (Wang 
& Yang  2001). 

Seven  sub-ranges  are  identified  with  each  having  distinct  dynamic  re- 
sponses to  the  finite  random  disturbances.  The  first  ranges  from  L2  = — 20 
to  L2  = —14.5,  where  the  finite  random  disturbances  lead  all  steady  solu- 
tions at  any  fixed  L2  to  a steady  symmetric  multi-cell  state  on  ASla  with 
the  same  L2.  The  second  covers  the  range  —14.5  < Dk  < —13.6  where 
there  is  co-existence  of  three  stable  steady  symmetric  multi-cell  states.  In 
the  third  sub-range  —13.6  < L2  < 12.1,  all  steady  solutions  evolve  to 
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a temporal  periodic  solution.  The  fourth  sub-range  is  from  L2  = —12.1 
to  L2  — —11.5  where  the  solutions  response  to  the  finite  random  dis- 
turbances in  the  form  of  either  periodic  oscillation  or  chaotic  oscillation. 
There  is  the  co-existence  of  periodic  and  chaotic  oscillations.  In  the  fifth 
sub-range  —11.5  < L2  < —10.5,  all  steady  solutions  evolve  to  a tempo- 
ral chaotic  solution.  The  next  sub-range  —10.5  < L2  < —10.2  serves  as 
a transition  between  the  chaotic  oscillation  and  the  stable  steady  2-cell 
flow.  The  solutions  response  to  the  finite  random  disturbances  in  the  form 
of  subharmonic-bifurcation-driven  asymmetric  oscillation.  In  the  last  sub- 
range L2  > —10.2,  the  finite  random  disturbances  lead  all  steady  solutions 
at  any  fixed  L2  to  a stable  steady  symmetric  2-cell  state  on  ASls  with  the 
same  L2.  A detailed  discussion  of  stability  features  can  be  found  in  Yang 
(2001). 


4.  Concluding  Remarks 

The  governing  differential  equations  from  the  conservation  laws  are  dis- 
cretized by  the  finite  volume  method  to  obtain  discretization  equations, 
a set  of  nonlinear  algebraic  equations.  The  discretization  equations  are 
solved  for  parameter-dependence  of  flow  and  temperature  fields  by  the 
Euler-Newton  continuation  with  the  solution  branches  parameterized  by 
L2,  the  arclength  or  the  local  variable.  The  bifurcation  points  are  detected 
by  the  test  function.  The  branch  switching  is  made  by  a scheme  approxi- 
mating the  difference  between  branches  proposed.  One  symmetric  and  two 
symmetric/asymmetric  solution  branches  are  found  with  fourteen  bifurca- 
tion and  seventy-five  limit  points.  Both  solution  and  flow  structures  are 
much  more  richer  than  those  available  in  the  literature. 

The  dynamic  responses  of  multiple  solutions  to  the  2D  finite  random 
disturbances  are  examined  by  the  direct  transient  computation.  The  finite 
random  disturbances  are  found  to  lead  the  steady  solutions  to  a stable 
steady  multi-cell  state  in  —20  < Dk  < —14.5,  the  co-existence  of  three 
stable  steady  multi-cell  states  in  —14.5  < L2  < —13.6,  a temporal  periodic 
oscillation  in  —13.6  < Dk  < —12.1,  the  co-existence  of  periodic  and  chaotic 
oscillating  states  in  —12.1  < L2  < —11.5,  a chaotic  oscillation  in  —11.5  < 
L2  < —10.5,  a subharmonic-bifurcation-driven  asymmetric  oscillation  in 
— 10.5  < L2  < —10.2,  and  a stable  steady  2-cell  state  in  —10.2  < L2  < 70. 
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Figure  1. 


Physical  problem  and  coordinate  system 
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(c)  -12.1  <L2<  -11.4 


(d)  -11.5  <L2<  -11.4 


Figure  2.  Solution  branches  and  limit/bifurcation  points  (a  = 0.02,  Pr=0.7,  Dk  = 300 
and  LI  = 28) 
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Figure  3.  Flow  and  temperature  fields  (a  = 0.02,  Pr=0.7,  Dk  = 300  and  LI  = 28;  left: 
secondary  flow,  middle:  streamwise  velocity,  right:  temperature) 


Figure  4-  Typical  secondary  flow  patterns  of  six  solutions  among  thirty-nine  solutions 
at  L2  = —11.7  (<7  = 0.02,  Pr=0.7,  Dk  = 300  and  LI  = 28) 


